Commuting Heisenberg operators as the quantum response problem: 
Time-normal averages in the truncated Wigner representation. 
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The applicability of the so-called truncated Wigner approximation (-W) is extended to multitime 
averages of Heisenberg field operators. This task splits naturally in two. Firstly, what class of 
multitime averages the -W approximates, and, secondly, how to proceed if the average in question 
does not belong to this class. To answer the first question we develop an (in principle, exact) 
path-integral approach in phase-space based on the symmetric (Weyl) ordering of creation and 
annihilation operators. These techniques calculate a new class of averages which we call time- 
symmetric. The -W equations emerge as an approximation within this path-integral techniques. 
We then show that the answer to the second question is associated with response properties of the 
system. In fact, for two-time averages Kubo's renowned formula relating the linear response function 
to two-time commutators suffices. The -W is trivially generalised to the response properties of the 
system allowing one to calculate approximate time-normally ordered two-time correlation functions 
with surprising ease. The techniques we develop are demonstrated for the Bose-Hubbard model. 

PACS numbers: 02.50.Ey,03.65.Sq,05.10.Gg 



I. INTRODUCTION 



One of the most fascinating areas of ultra-cold atomic 
physics is the experimental investigation of the dynam- 
ical properties of interacting many-body systems. The 
control of experimental parameters and ability to tailor 
systems is allowing many interesting effects to be ob- 
served which would have been almost impossible in the 
recent past. Among important examples are the recently 
observed dynamical instabilities [l[ and inhibited trans- 
port m one dimensional (ID) and three-dimensional 
(3D) optical lattice systems, as well as nonlinear self- 
trapping in ID periodic potentials [ij and in Josephson 
junctions Despite these important experimental ad- 
vances, the theoretical description of dynamical proper- 
ties, especially for strong interactions, remains a major 
challenge, with progress having been incremental up to 
now. For bosons, except for rare cases where the dynam- 
ics are analytically tractable, e.g. by Bethe ansatz Q, 
one method is to adapt the phase-space representations of 
quantum optics @) IS & B3 ■ Attempts have been made, 
for example, to apply the positive-P representation to the 
dynamics of trapped and colliding Bose-Einstein conden- 
sates [111 El EH [lj] , although these have been only par- 
tially successful due to numerical instabilities which mean 
that integration is restricted to either short times or small 
interaction strengths. We note here that efforts are being 
made to extend the usefulness of this method [HI, [H, [TlJ , 
with promising results for single-mode systems. For ID 
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systems numerical algorithms based on adaptive matrix- 
product decompositions of the state vector, such as the 
t-DMRG [H[ and the related TEBD algorithm have 
recently been developed. They too are restricted how- 
ever to short times or systems with slow entanglement 
growth. 

For highly nonlinear undcrdamped systems such as 
trapped condensates, it is often easier to use the trun- 
cated Wigner approximation (-W) [20|] , with the main ad- 
vantage being that it is numerically stable [3, HH HH 
and simple to implement. The -W takes into account 
initial uncertainty between conjugate variables or initial 
quantum noise [24j . which is often necessary to trigger 
certain dynamical processes. This method was recently 
used to predict and explain several experiments in the 
context of interacting cold atom systems. As an exam- 
ple, using this method the damping of a dipolar mo- 
tion of a ID condensate in an optical lattice was first 
predicted in ref. [25[ and later results were qualitatively 
confirmed in experiment 0. This experiment was later 
simulated more accurately using a multi-band version of 
the -W f26{. In Refs. [23, [n[ the -W was used to ana- 
lyze splitting or merging between elongated condensates, 
closely mimicking the situation realized in recent exper- 
iments 29]. In Ref. 30], the -W enabled the authors to 
explain the coherence dynamics after a sudden quench of 
tunneling from an insulator to a superfluid, giving good 
agreement with experimental results [3l| . For a fuller re- 
view of other recent developments, we direct the reader 
to Ref. H3. 

While the truncated Wigner representation is becom- 
ing increasingly utilized, there are two drawbacks which 
prevent it from being the numerical technique of choice. 
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The first is its approximate nature, arising from the trun- 
cation procedure, which may sometimes lead to demon- 
strably wrong results [33, 34]. In principle, this drawback 
can be overcome. A way to fully map the quantum mas- 
ter equation onto Wigner representation stochastic differ- 
ence equations has been developed, but does not result 
in a widely and easily applicable method [H, HH, H3] • An 
expansion allowing one to include the dropped contribu- 
tion perturbatively via quantum jumps was suggested in 
Ref. (24|, but until now this has been restricted to the 
calculation of single-time averages. 

The second drawback is that averages of the phase- 
space variables in the Wigner representation do not map 
directly onto expectation values of time-normally ordered 
multitime operator products corresponding with experi- 
mental measurements. For two operators, we recover a 
symmetrised product, while for three and more opera- 
tors, there emerges a new type of ordering of Heisenberg 
operators which we term time- symmetric [361 ]. Convert- 
ing these to time-normal ordering requires the calcula- 
tion of commutators of Heisenberg operators at differ- 
ent times. At first glance this drawback appears to be 
fundamental, but we will demonstrate that it can be cir- 
cumvented with surprising ease, by making use of the 
fact that commutators of Heisenberg operators at dif- 
ferent times express the response properties of a quan- 
tum system [11, [M HI Efl EH, Si- Generally, com- 
muting the Heisenberg operators requires solving the full 
quantum stochastic nonlinear response problem, and, al- 
though possible in principle, this remains a formidable 
task [42l. [43(] . For two-time commutators, however, the 
process is simplified due to Kubo's linear response formu- 
lation [38l [3^ | , which has previously been used to convert 
normally-ordered correlations into symmetrically ordered 
ones within the positive-P representation [34|, |44| . In 
this paper we show how, combined with the truncated 
Wigner representation as a computational tool, Kubo's 
linear response theory turns into a simple yet powerful 
approximate method of calculating two-time correlation 
functions of interacting bosonic fields with arbitrary op- 
erator ordering. We stress here that the Kubo relation 
is exact, with the approximate nature coming from the 
truncation process used to find the appropriate stochastic 
equations. 

In more traditional phase-space techniques, the 
truncated- Wigner equations are developed by dropping 
the third-order derivatives in the generalised Fokker- 
Planck equation for the single-time Wigner distribution 
floL [20( 1 . The corresponding Langevin equations are ei- 
ther non-stochastic (without losses) or probabilistic (with 
losses) and are very easy to handle numerically. However 
simple and straightforward, this conventional way of de- 
riving the truncated Wigner approach leaves it unclear 
if it can be applied to any multitime quantum averages. 
In our approach, the truncated- Wigner equations emerge 
as an approximation within rigorous phase-space path- 
integral techniques. By itself, the path integral expresses 
the averages of time-symmetrically-ordered products of 



Heisenberg operators. The generalisation associated with 
extending the truncated- Wigner equations to multitime 
averages is thus highly nontrivial and requires a new con- 
cept: the time-symmetric ordering of the Heisenberg op- 
erators. There does not seem to be a way of guessing this 
concept from within the conventional phase-space tech- 
niques. Just proving the equivalence between the time- 
symmetric and the conventional symmetric ordering of 
free-field operators is a nontrivial task [43| . 

To make this paper more accessible we will begin with 
a theoretical summary in section [TU where we list the 
key results, supporting them by leading considerations. 
The actual theory is presented in section IIIIl Sections 
IIVI and fVl illustrate how to apply the method in practice. 
The summary of section [H] suffices for understanding the 
examples in sections IIVI and [Vj The reader who is not 
interested in a rigorous justification of the method can 
safely ignore it. Our first example is Kerr oscillator (sec- 
tion |IV|). This is an exactly soluble problem; moreover, 
all calculations in the truncated Wigner approach can 
be carried out analytically. In section [V] we apply the 
method to the Bose-Hubbard model, comparing the re- 
sults to direct calculations in Hilbert space. We consider 
relatively small chains consisting of few sites to make ex- 
act calculations possible, although the truncated Wigner 
approach can be extended to much larger systems (see 
e.g. Refs. [ill Hfjj]). We find good agreement between 
exact and approximate results over relatively short time 
scales, with the applicability of this method to longer 
times remaining an open problem. 



II. THEORETICAL SUMMARY 

A. Symmetric representation of Heisenberg 
operators 

We will assume that the reader is familiar with the 
phase-space basics, including the concepts of symmet- 
ric (Weyl's) operator ordering, symmetric representation 
of operators and the Wigner function. The necessary 
minimum of information is summarised in section IIII Bl 
Details may be found in Refs. [H H H H| . 

The formal techniques we develop in this paper extend 
the well-known symmetric representation of quantum op- 
tics [Io| to multitime problems. First and foremost, in 
place of a quasiaverage over the Wigner quasiprobability 
distribution we find a phase-space path integral. Within 
the path- integral techniques, we derive generalised phase- 
space correspondenses mapping multiplication of a q- 
number quantity by a creation or annihilation operator to 
phase-space. Conventional phase-space correspondenses 
apply to free operators and the generalised ones apply 
to Heisenberg operators. They are exact and do not de- 
pend on the nature of nonlinearity (interaction). Simi- 
larly to the way in which conventional phase-space cor- 
respondences allow one to reorder creation and annihila- 
tion operators, the generalised ones allow one to reorder 
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a pair (say) of Heisenberg operators with unequal time 
arguments. These techniques are not restricted to any 
special Hamiltonian and can be employed for all bosonic 
systems. 

To be specific, consider an illustrative example of the 
anharmonic oscillator (Kerr oscillator) with the Hamil- 
tonian 



~ Tin ^+2-2 
Hq — — a< a . 
2 



(1) 



Here, a,a^ is the standard creation/annihilation pair, 

[a,ot]=i. (2) 

The Heisenberg "field operators" are defined in the nor- 
mal manner as 



A(t) =U\t,t )aU(t,t ), 

A\t) =wt(i,t )attf(t,to) 

Here, U(t,to) is the evolution operator, 
LA [t, to) — exp 



(3) 



(4) 



where to is the coincidence point for the Schrodinger and 
Heisenberg pictures. For an arbitrary Ho the evolution 
operator is introduced through the Schrodinger equation, 



ifi 



dU(t,t a ) 
dt 



H U(t,t ), U(t ,t ) = l. 



(5) 



Note that we do not divide the Hamiltonian into the 
free and interaction part, nor introduce the interaction 
picture, nor free-field operators. The methods we develop 
in this paper are strictly nonperturbative. The truncated 
Wigner representation does not correspond to linearisa- 
tion of any kind. It follows separation of nonlinearity 
and noise for path-integral trajectories with neglecting 
the latter, for details see section IIII Dl There does not 
seem to be a way of even formulating such approxima- 
tions in Hilbert-space terms. 

In section Hill we construct a phase-space path- integral 
approach which allows one to calculate the quantum av- 
erage of the symmetrised product 



{T w A\h)A(tx)) = a(ti)tt*(t 2 ), 



where 



T w A\t 2 )A{t x ) = ^\A\t 2 )A{t l ) + A{h)A^t 2 ) 



(6) 



(7) 



For the time being, Tw is just an ad hoc notation; its 
actual meaning will be the subject of section III CI The 
quantum averaging is over the Heisenberg p-matrix p, 
(with X being an arbitrary operator) 



(X) = TrpX. 



(8) 



The double bar denotes symbolically the path integral 
over the c-number trajectories a(t). It can be thought 
of as a stochastic average with a nonpositive "measure" , 
whose exact meaning will be clarified in the next section. 
What is important is that the path integral ([6]) is ap- 
proximated by the truncated Wigner approach. In this 
approximation the trajectories a(t) are deterministic and 
obey the equation 



,da(t) 
1 dt 



i(t)\ 2 - l]a(i). 



(9) 



This equation follows both from the traditional phase- 
space methods [20( and from the path-integral tech- 
niques in section IIIIl The trajectories being determin- 
istic, stochasticity in Eq. ( reduces to the averaging 
over the Wigner function W(a) corresponding to the 
Heisenberg p-matrix (assuming this function is nonncg- 
ative, W(a) > 0). The path integral in Eq. ( [6|) then 
reduces to averaging over solutions of Eq. ( [5]) with a 
random initial condition, 



a(ti)a*(t 2 ) - / —W(a)a(ti)a*(jh), (10) 

where a(t) is a solution to (|15p with the initial condition 
&(to) = ct. Making trajectories deterministic is an ap- 
proximation, so that the average (|10p only approximates 
the quantum average, 



(T w A\h)A(h)) w a{tx)a*{t 2 ), 



(11) 



unlike the path integral © which is exact. 

It is worth stressing here that the results of this paper 
are not Eqs. fl9j)- (fTTj) as such, but the fact that they ap- 
ply with t\ 7^ t% . The truncated Wigner approach was de- 
rived within the conventional phase-space techniques. By 
construction, it is only applicable to time-dependent av- 
erages of Schrodinger operators, or, equivalently, to aver- 
ages of Heisenberg operators with equal time arguments. 
In other words, conventional phase-space methods allow 
one to verify Eqs. ( l5])- (rTTj) only for t\ = t 2 . Exten- 
sion to unequal times (physically, to spectral properties 
of the system) requires alternative techniques such as a 
phase-space path integral. 



B. Generalised phase-space correspondences 

Our next goal is to extend the truncated- Wigner fur- 
ther by lifting the ordering restriction of Eq. ( [5]). 
Namely, we wish to calculate the time-normally ordered 
average, 

(AHh)A(t 2 )) = (T w Al(t 2 )A(ti)) 

-i<[i(ta), .^(ti)]), (12) 

where it has been expressed by the symmetrised average 
and the commutator. The symmetrised term is given by 
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([6]). To express the commutator we employ the same 
idea as was used in Ref. [44j to reorder a time-normal 
average symmetrically. Namely, we assume that t 2 > t-y 
and relate the commutator to the linear response of the 
system, 



S(A(t 2 )) 



Ss(h) 



(13) 



This relation is simply Kubo's formula for the linear re- 
sponse function [38|, |39( written "from right to left." It 
implies that the Hamiltonian of the system has been com- 
plemented by an interaction with the external c-number 
source s(t), 



Ho -> H(t) =H Q - s{t)a) - s*(t)a. 



(14) 



Strictly speaking, the condition s — must then be ap- 
plied to both sides of Eq. ( I13p (and in fact to all quantum 
averages in the above), but in practice it suffices to re- 
member that the only quantity defined with s =/= is 
(A(t 2 )) in Eq. (E3J). 

For simplicity we will confine the rest of the discus- 
sion to the truncated Wigner representation. In this case 
the trajectories a(t) are deterministic and obey the equa- 
tion which differs from §§§ by the presence of an additive 
source: 



da(t) 
dt 



n[\a{t)\ 2 - l]a(f) - s(t). 



(15) 



Within conventional phase-space methods, we find this 
equation by noting that the linear interaction terms in 
the Hamiltonian contribute only to drift terms in the gen- 



eralised Fokker-Planck equation. A path-integral deriva- 
tion of (|15p extending its applicability to multitime av- 
erages will be given in section MI Dl 

To calculate the linear response function (|14[) one needs 
an infinitesimal instantaneous source at t = t\, 



s(t) = —ihSaS(t — ti). 



(16) 



Substituting this source into Eq. ( [15|) we see that it 
causes a discontinuity of the trajectory at t = ti (a 
"quantum jump" in the terminology of Ref. [24]). The 
additional factors in (116|) were chosen so as to make this 
discontinuity exactly equal to 8a. We thus have a simple 
correspondence between sources and "quantum jumps" : 



d 



d 



Ss(ti) hda(ti)' 6s*(ti) 



Kda*(h) 



(17) 



Mathematically, the derivatives d/da(ti), d/da*(t\) cor- 
respond to a variation of the initial condition set at t = t\ 
instead of t — to. Such notation is to some extent infor- 
mal but convenient. For the commutator we then find 
(*i < h) 



([A(t 2 ),A\h)]) 



da(t 2 ) 
da(t\) 



t 2 > ti. 



(18) 



This relation follows from the correspondences (jT7|) and 
from the path-integral representation of the average. 



{A{h)) = a(t 2 ). 



(19) 



Using (|18j) we can express the time-normal average in the 
Wigner representation, 



!a*(ti)a(t 2 ] 

a*(t!)a(t 2 ) 

r 



The second line here follows by conjugating the first one 
and replacing t\ <-> t 2 . 

These generalised phase-space correspondences (|20p are 
the central result of the paper. Certainly, the above 
derivation is no more than leading considerations, but 
the rigorous treatment in section ITlI Fl gives the same re- 
sult. Moreover, it shows that Eqs. ( |2U|) hold not only 
as an approximation in the truncated Wigner represen- 
tation, but also as an exact relation within the rigorous 
path-integral approach — in which case the bar in ([2TJ]) 
should be replaced by double bar. 

Why do we call Eqs. ( |2"0|) "generalised phase-space 
correspondences?" To recognise the connection we take 
the limit t±,t 2 — > to- Using the fact that that A(to) = 



a, A^(to) — d^ and dropping the time argument in a(to) 
we find 




where the averaging is simply over the Wigner function 
W(a), cf. Eq. ([10]). First of all, both relations in Eq. ( 
I2T1) are correct. Indeed, they result in the formula 

(tfa)=W-l- ( 22 ) 
The average over the Wigner function expresses symmet- 
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rically ordered products of a, at; in particular, 

j^pr = I (aat + a t a) = (at a) + 1, (23) 

in obvious agreement with (|22[) . Furthermore, Eqs. ( |2"T|) 
are particular cases of the standard phase-space corre- 
spondences, 



H - (H^) r <°>- <24) 
< a,f ->=(«*-^) F<£,) - (25 » 

where 1" is an operator and Y(a) is its symmetric repre- 
sentation. The first of Eqs. ( 121]) follows from Eq. (]2JD 
with r = a^, F(a) = a* , and the second one — from Eq. 
( [23]) with y = a, y(a) = a. Multitime generalisations 
of Eqs. (|2"3|). ([25]) are derived in section ImFl 

We conclude this paragraph with a remark on termi- 
nology. To maintain rigor, one should distinguish shifts 
of trajectories effected by instantaneous sources (|16l) from 
quantum jumps. The latter term was introduced in Ref. 
[2J|, where discontinuities of trajectories were used as 
formal means to express perturbative corrections to the 
truncated Wigner approach. These corrections come 
from quantum noises which do not have any classical 
interpretation whatsoever, while the c-number external 
source is to a large extent a classical object. Maintaining 
the distinction between shifts and quantum jumps thus 
appears physically justified. However, such clear-cut dis- 
tinction is an artifact of an undamped model with quar- 
tic interaction. For instance, for the damped harmonic 
oscillator the equation for the Wigner function is a gen- 
uine Fokker-Planck equation. In this case the "quantum 
noise" is fully probabilistic, i.e., classical. We will use 
"quantum jump" as a blanket term applicable to both 
types of discontinuities. 



C. The time-symmetric ordering 

The fact that the path integral © calculates (and the 
truncated Wigner approach approximates) symmetrised 
products of Heiscnberg operators does not generalise 
to products of three and more operators. Instead of 
fully symmetrised products, one discovers a new type 
of ordering of Heisenberg operators, which we call time- 
symmetric and denote as Tw- We find this interesting and 
important enough to be worth reporting, notwithstand- 
ing the fact that it is not directly relevant to purposes of 
this paper. 

A time-symmetrically ordered product of the "field op- 



erators" Ait) , A* (t) is defined recursively as 
T w i = i, r w A{t) = A{t), T w AHt) = A\t), 
T w A{t)V [>t] = \{A{t),T w V [>t] ), (26) 
T w A\t)T [>t] = \{A\t),r w T\ >t] }. 

Here, 'P[>t] is a product of field operators with all time 
arguments exceeding t; the curly brackets stand for the 
anticommutator, {X, y} = Xy + yX. It is implied that 
under the sign of "TV-ordering the field operators com- 
mute freely. The quantum average of an arbitrary time- 
symmetric product is expressed as a path-integral aver- 
age, (m, n > 0) 

(T w A(h) ■■ -A{t m )A («£)••■ it (4)) 

= a(t 1 )---a(t m )a*{t[)---a*(t' n ). (27) 

For the exact meaning of this relation we refer the reader 
to the section Hill Again, what matters is that the trun- 
cated Wigner approach represents the path-integral av- 
erage approximately, 

{T w A{t l )---A{t m )AHt' 1 )---AHt' n )) 

« a(h) ■ ■ ■ a(t m )a* (*!)••• a* (t' n ) ^ 

= / ^W(a)a(t 1 )---a(t m ) a *(t' 1 )---a*(t' n ), 

cf. Eq. (HP]). Equations ([26]) and ([27]) may be directly 
generalised to multimode and real-space cases, by supple- 
menting the time arguments by suitable "labels," such as 
mode indices or spatial arguments. For an example (the 
Bose-Hubbard chain) see section [V] 

The two most important properties of the time- 
symmetric products are: these products are continuous 
at coinciding time arguments, and for free-field opera- 
tors they turn into the conventional symmetric (Weyl) 
ordered products. For two operators, the time-symmetric 
product coincides with a symmetrised product given by 
Eq. ([7]). That quantity is continuous at t — t'\ more- 
over, for coinciding times, ([7]) reduces to the conventional 
formula for the symmetric ordering, which naturally ap- 
pears in the -W approximation [32, IH, [49j] , 

Wjaat} = i(aa f +a)a). (29) 

(Recall that for coinciding times the field operators com- 
mute the same way as the creation and annihilation op- 
erators.) For three operators and t\ < t% < t% we have, 
for example, 

TwA{ti)A{t 2 )A\h) 

= -^[AitM^AHh) + A{t 2 )AHh)A(h) 
+ MtMHt 3 )Mt2) + i t (i 3 )i(*2)i(ti)] . (30) 
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Here the time-symmetric ordering is not the same as the 
fully symmetric ordering; in the latter there should be 
two additional terms with A{t\) in the middle. Again, it 
may be shown that (|30[) is continuous at coinciding time 
arguments, and that with all three times equal it agrees 
with the formula for the Weyl-ordered product, 

W{a 2 a f } = - (a 2 a t + a t a 2 + aa t a) . (31) 

o 

Detailed discussion of the time-symmetric ordering re- 
quires advanced formal tools and will be presented else- 
where [43| . 

We note that all operator products entering the time- 
symmetric product exhibit a special order of time argu- 
ments: times first increase then decrease. Such order 
of operators is characteristic of Schwinger's closed-time- 
loop formalism [5l|. This connection is investigated in 
Ref. [43j]. We also note without proof that only such 
"Schwinger-ordered" operator products have causal rep- 
resentation through quantum jumps similar to Eqs. ( I20p . 
This restriction becomes nontrivial for products of three 
or more operators. For example there is no causal repre- 
sentation through the response for finding the expecta- 
tion value of A(t2)A(ti)A(t 3 ) with t\ < t 2 ,t 3 . For this 
particular ordering one cannot avoid finding the response 
at t\ to a perturbation which happens later in the evo- 
lution either at t = t 2 or t = t 3 . For more details see 
Ref. \E^. 



III. MULTITIME WIGNER APPROACH 
A. Preliminary remarks 

In this section we present a rigorous derivation of the 
"generalised phase-space correspondences" (|20|) . The 
reader who is interested only in applications of the 
method can safely skip the formalism and go directly to 
examples in sections IIVI and [V] 

For simplicity we will continue working with the il- 
lustrative example of the Kerr oscillator. The necessary 
definitions were given in section Hi Al In fact all formulae 
in this section apply to arbitrary time-dependent Hamil- 
tonians, and can also be easily generalised to multimodc 
problems, simply by complementing the time arguments 
by other "labels," such as mode indices or spatial argu- 
ments. 



B. Phase-space basics 

For the reader's convenience, we summarise here the 
necessary facts from phase-space techniques [13, IH, |49|, 
[Hoj . The displacement operator is defined as 

D{a) =e a ^'- a * d . (32) 



For an arbitrary operator A, one introduces its charac- 
teristic function, 

XA (a) = Tri£)t(a), (33) 
and its symmetric representation, 

A(a) = [A] (a) = J f£ XA ((3)e^'^' a . (34) 
Expressions for A in terms of these read 

I , (35) 

The notation [• • • ] (a) is convenient for symmetric rep- 
resentations of operator expressions, as, for instance, in 
Eqs. ([35]), (gOJ) below (see also endnote HH). Of use to 
us will be the relations, 

[a] (a) = a, [a^] (a) = a* , 

[a t a](a) = |a| 2 -i, (36) 

[a t2 a 2 ](a) = \a\ 4 - 2\a\ 2 + ^. 

Displacement operators form a complete set with re- 
spect to the Hilbert-Schmidt norm, 

/,2 
— TrAD(a)TTD\a)B 

= / —A(a)B(a). (37) 

The last equation here is a consequence of (|34j) . In partic- 
ular, it allows one to write a phase-space representation 
of a quantum average, 

/i'2 
—A{a)p{a). (38) 

Of importance to us will be a relation expressing the 
Wigner representation of an operator product AB by the 
Wigner representations of the factors: 

[AB](a)= J rf2 ^% («-aoK-(<*-<*o)V 

x A(a )B(a + a/2). (39) 

By the change of variable ao — > ao + a/2 we can write 
this in the alternative form, 

[AB](a)=f 

x A{a Q - a/2)B{a Q ). (40) 
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These relations follow from expressing the operators by 
their symmetric representations using Eq. ( I35[) and then 
employing Eq. ( [34j) to express L473](a). It is easy to 
verify that [501 ] 

TrD(a)D((3)D(/3') 

= ir5W(a + p + f3')e^"-W'\ (41) 

The rest of the calculation leading to Eqs. ( [35]), ([4*0]) is 
straightforward. 

The symmetric representation of an operator is often 
introduced as an expression for this operator in terms of 
symmetrically (Weyl) ordered products of creation and 
annihilation operators. Such products are defined postu- 
lating that 

[W{a m a tn }] (a) = a m a* n . (42) 
Equations (|36[) are then written as operator formulae, 
a = W{a}, a^Wja 1 }, 

(43) 



a' a 



W{aot}--, 



at 2 a 2 = W{a 2 at2}_2W{aat} + 1 



These equations (and thus Eqs. ( 136"]) ) may be verified 
noting that the displacement operator is naturally Weyl- 
ordered, 



D(a) =W{D(a)}, 



(44) 



and can therefore serve as an operator-valued character- 
istic function for the symmetrically ordered products, 



D(a) 



y an{ ~?, )m w{a m tf n }. 

±- — J 771 ' 11 ' 



(45) 



Verification of Eqs. ( [43]) reduces to developing D(a) in 
a power series, with the subsequent use of {2]). 



C. Phase-space transition amplitude 



With the only exception of Eq. ( 132]) which employs the 
creation/annihilation pair in the Schrodinger picture, all 
definitions in scction llll Bl mav be applied to Schrodinger 
as well as to Heisenberg operators. If a particular op- 
erator is time-dependent, its symmetric representation 
is also time-dependent. The time-dependent symmetric 
representation of a Schrodinger operator and the time- 
dependent Wigner function are defined as follows, 

[73(f)] (a) = B(a, f), [p(t)] («) = p(a, t), (46) 

cf. Eq. ( 134"]) . We stress that both definitions here are for 
operators in the Schrodinger picture. In the Heisenberg 
picture the density matrix is stationary and coincides 



with p(t ); its symmetric representation thus coincides 
with p(a, to)- The Heisenberg counterpart of 73(f) reads 



23(f) = W(t,t )B{t)U(t,t ). 



(47) 



We do not introduce any special notation for symmet- 
ric representations of Heisenberg operators but use the 
bracket symbol instead, cf. Eq. ( |5"T]) below. 

Using Eq. ([37]), the quantum average of 23(f) may be 
written as 



(23(f)) = / ^B(a,t)[U(t,t )p(t )W(t,t )](a) 
cPctocPa 

= — B{a, t)U(a, f, a ,to)p(a , f ), 



(48) 



where we have introduced the phase-space transition am- 
plitude 

U(a,t,a ,t ) = J *P°fP e <*(r-«'fi+aoK-eSfh 

x TTD\f3)U(t,t )D\p )U^t,t Q ). (49) 

By construction, this amplitude evolves the Wigner func- 
tion in time, 



d 2 a 



U(a,t,a ,t )p(a ,t ) , (50) 



but it can also be applied to the operator, 
f d 2 

0t o (t)](a o )= I -^B(a,t)U(a,t,a ,t ). (51) 

In this formula the dependence of the Heisenberg opera- 
tor on the coincidence point fo is made explicit showing 
it as a subscript. Such notation is convenient when the 
coincidence point itself becomes a variable as in section 

[The] below. 



D. Phase-space path integral and the truncated 
Wigner representation 

The group property of the evolution operator 

U(t,t ) = L>(f,fi)L>(fi,f ), f>fi>f , (52) 

results in the related property of the transition ampli- 
tude, 



U(a,t,a ,t ) 



d 2 ai 



U(a,t,a 1 ,ti)U{a 1 ,ti,a ,to) ■ 



(53) 



Breaking the time interval [fo, f] into A7 + 1 Trotter slices, 
f-f 



Af 



M + 1 



, f fe = f + A;Af, k = 0,...,M, (54) 
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we can define the path-integral representation of the 
phase-space amplitude as the limit 



U(a,t, a ,t ) = lim U(a,t,a M ,tM) 

M—*oo 



n 

k=l 



d 2 a k 



U(ak,tk,ak-i,tk-i) ■ (55) 



Each amplitude on the RHS here is over an infinitesimal 
time interval At. 

To understand the path integral we have thus to un- 
derstand the infinitesimal transition amplitude. It may 
be evaluated using the method introduced by one of us 
in Ref. [24| . We start from the von-Neuman equation for 
the density matrix 



ihfo) = [H(t),p(t)l 



(56) 



so that 



iAt - iAt 
p{t + At) = p(t) - —H{t)p{t) + —p(t)H(t). (57) 



Note that we wrote H (t) to highlight the fact that the 
derivation is valid for arbitrary time-dependent Hamilto- 
nians, including Hamiltonians with external sources such 
as (flip. Employing Eqs. ( El]), (SOI) and introducing 
the symmetric representation of the Hamiltonian in the 
Schrodinger picture, 



H(a,t) = [H(t)}(a), 



(58) 



we have 



p(a,t + At) 



7T- 



fr(ao-|,i) -fr(a + |,* 



(59) 



see also endnote (5|| . Comparing this to Eq. ( [50]) and using the fact that At is infinitesimally small we find 
U (a, t + At, ao, t) 



d 2 a 



■ exp 



At 

a- a +if(a ,t) — 



At' 

a - ao + if(ao,t) — 



t iAt 

a+ l —h^(a ,a,t)\, (60) 



where f(ao,t) and (ao, cr, t) are found by expanding 
the symmetric representation of the interaction Hamilto- 
nian into power series: 



f(uo,t) 



dH{a ,t) 



f*(a ,t) 



dH(a ,t) 



dao 



h ( - 3 \a ,a,t) = H^a 

- o-f*(a ,t) - a*f(a ,t) 



lA-Hlao-lt 



(61) 



The term (ao , cr, t) is responsible for cubic noise, 
which accounts for quantum fluctuations; a consistent 
derivation of the path integral with the cubic noise will 
be subject of a separate paper. Attempts to simulate 
the cubic noise numerically were rather disappointing 
[35l HH, H3- I n R- er - 01 one 01 us showed how it can 
be taken into account perturbatively through the non- 
linear response. In Ref. [46l ] this nonlinear response was 
implemented to improve the accuracy of the -W for a 
large BH chain of 128 sites. 

On the other hand, neglecting cubic noises simplifies 
our task enormously by removing all mathematical prob- 
lems associated with their highly singular nature. With- 
out the integral in (|6"0|) is calculated straightaway, 



and we have 

U(a,t + At, a a , t ) = ir5 {2) ( a — ao + Trf{oto, t)At 

(62) 

This corresponds to a deterministic evolution in phase- 
space along the trajectories satisfying the equation 



iha = f{a, t) . 



(63) 



By making use of Eqs. ( |3"6"|) . for the Kerr oscillator 
we find Eq. ( [9]). We have thus recovered the well- 
known truncated Wigner representation (20| . However, 
unlike in Ref. [2(| , we have found it as an approximation 
within a consistent phase-space path-integral approach. 
This allows us to answer two questions which cannot be 
answered in the derivation based on the Fokker-Planck 
equation. Firstly, which quantum averages the path in- 
tegral calculates, and, secondly, how one could evaluate 
other types of averages. This will be subject of Secs. lIII El 
andniFl 

Equations (|6ip make it obvious that external sources in 
the Hamiltonian manifest themselves as additive sources 
in the equations for trajectories. Indeed, using Eqs. ( 
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|36|) . for the Hamitonian ([14"]) we have, 

[H(t)] (a) = [Ho] (a) ~ s{t)a* - s*(t)a. (64) 
The source terms only modify the regular evolution, 



/(<*,«)-> f{a,t)-s(t), 
/*(a,t) ->/*(«,*) -«*(«), 



(65) 



For the Kerr oscillator this results in Eqs. ([15]). Equation 
([64]) holds for an arbitrary iSo, so that replacements ([65]) 
apply in general. 

That the Kubo-style sources in the Hamiltonian ap- 
pear as additive sources in the equations of motion for 
the phase-space trajectories is in fact true for arbitrary 
phase-space techniques. Indeed, irrespective of the oper- 
ator ordering, linear terms in the Hamiltonian manifest 



themselves only as drift terms in the generalised Fokker- 
Planck equation and thus only as additive terms in the 
corresponding generalised Langevin equations. For an 
example see Ref . [44| , where external sources were intro- 
duced in the positive-P representation. 



E. Time-symmetric operator ordering 



We will now address the question of which quan- 
tum averages the path integral calculates. To define 
this more clearly, consider the path-integral average 
(ti < h < ■ ■ ■ < t K ) 



a(ti)a(t 2 ) ■ --a^K) 



cPaocPai ■ ■ ■ cPctK 



x a K U(aK,tK,a K -i,t K -i)aK-iU(a K -i,t K -i,a K -2,t K -2) ■ • • aiU(ai, ti, a , t )p(a , t ). (66) 

We presume that there exists a rule of ordering for Heisenberg operators, which we term the time- symmetric ordering 
[36| and denote 7V, such that 



a{tx)a{t 2 ) ■ ■ ■ a(t K ) = (T w A(ti)A(t 2 ) ■ ■ ■ A{t K )). (67) 

The Heisenberg field operators are given by ([3]). It is easy to obtain a recursion relation expressing 
T w A{ti)A{t 2 ) ■ --Aitx) by T w A{t 2 ) ■ ■■A{t K )- Comparing Eqs. (155)1. (|57 )l to ([38 )1 we have 



[T w A to (ti)A to (<a) • • • A to (t K )] (a ) 



d 2 cti ■ ■ ■ d 2 otK 



x aKU(a K ,t K ,aK-i,tK-i)aK-iU(a K ^i,t K -i,a K -2,tK-2) • ■ ■ a\U{a\, t\, ao,*o), (68) 

see also endnote (5|| . In this relation the dependence of the Heisenberg operators on the coincidence point is made 
explicit. Applying it to the product TwAt 1 (t 2 ) ■ ■ ■ At^tx) with the coincidence point set at t\ we find 



[T w A tl (t 2 )--- A tl (t K )](ai) 



d 2 a 2 ■ ■ ■ d 2 otK 



x aKU(aK,tK,ceK-i,t K -i)aK-iU(a K ^i,t K ^i,aK-2,tK-2) • ■ ■ a 2 U(a 2 ,t 2 ,a 1 ,t 1 ). (69) 
Comparing Eqs. ( [55]) and ([M)> we see that 

f d 2 

[T w At (ti)A to (t2)---At (t K )](ao)= / a^ia^h, a ,t ) [TwA tl {t 2 ) ■ ■ -A^^k)] («i). (70) 
We now recall the standard phase-space correspondence, 

a[A](a) = l -[aA + Aa](a) = ±[{&,A}](a), (71) 
where the curly brackets stand for the anticommutator, {X, j 7 } = Xy + yX . This allows us to write 

ai [T w A tl (h) ■■■A tl (t K )] (ax) = ~ [{A tl (ti), T w A tl (t 2 ) ■ ■ ■ A tl (t K )}] (ax) , (72) 
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where we have used the fact that, with the coincidence point set at t = t\, the Heisenberg field operator A(t\) coincides 
with its Schrodingcr counterpart, 



A tl (ti) = a. 



(73) 



Equation (|70|) then becomes 

[T w A to {h)At {t 2 ) ■ ■■A t0 (t K )] (ao) 



1 f d 2 a x 



U(a u h, a , t ) [{A tl (ii),T w A tl (t 2 ) ■ ■ ■ A tl («*-)}] (<*i) ■ ( 74 ) 



2 J ir 

We now note that Eq. ( l5"Tj) is based solely on Eq. ( |4"7]) and is therefore a particular case of a more general relation 



[tii(t,t )Xti(t,to)](ao)= / —[X](a)U{a,t,ao,t Q ) 



(75) 



where the operator X may be arbitrary. Applying this to (|74|) we have 

[T w it (ti)i* (t 2 ) • ■•i i0 (i f f)](a ) = 2[W t (*i.*o){A 1 (ti),T' w A(*2) • ■ • K (t K )}U(ti,t )] (a Q ) 



; [{Ato (h),r w A t0 (t a ) • • • A (t* )}] (a ) . (76) 



We have thus arrived at the desired recursion relation, 



T w A(t 1 )A(t 2 )---A(t K ) = -{A(t 1 ),T w A{h)---A{t K )}, 



(77) 



where the dependence on the coincidence point has been 
dropped. 

If we replace a(ti) — > a*(ti) and A(t\) — > ^(ti) in 
([B"7]). equation ((771) will also hold with J.^) -> ^(ti). 
Furthermore, any subset of the factors a(t 2 ) • • ■ a{tx) 
may be complex-conjugated, provided the correspond- 
ing operators under the T^-ordering arc Hcrmitian- 
conjugated. As a result, we arrive at the recursive defi- 
nition of the time-symmetric ordering given by Eq, 



(175)) we have 

(y(t 2 )A(h)) = Try tl (t 2 )ap{h) = 
d 2 aid 2 a 2 



Y(a 2 )U{a 2 ,t 2 ,a 1 ,t 1 ) [ap(h)] (oi). (79) 



In this relation we made the dependence of y(t) on the 
in section QTC] For a brief discussion of this concept we coincidence point explicit, cf. Eqs. ( [70]), (J7JJ) and ([731) . 
refer the reader to that section. Detailed analyses will be The standard phase-space correspondences then allow us 



presented elsewhere [43 



to write 



F. Commuting Heisenberg operators as a response 
problem 



[ap(h)] (ai) - (ar + i^)p( Q!l , t t ). (80) 



Now we consider what happens if the quantum average 
we wish to calculate is not a time-symmetric one, but a 
time-normal one. In this paper, we only consider two- 
time averages (to < ti,t 2 ) 

{X{h)y(t 2 )) = Trpi^XihWih), (78) 

where X(t),y(t) = A(t),A'(t). A general discussion will 
be given elsewhere. Rather than distinguishing the cases 
t\ > t 2 and t\ < t 2 , we assume that t\ < t 2 and consider 
two distinct averages, {X(t\)y{t 2 )) and {y{t 2 )X{t\)). 
Consider, for example, the average 

<y(*a)^(*i)>. Moving 
the coincidence point to t = t% and using Eqs. ( |4"8|) and 



Using Eq. ([50]) to express p{a\,t\) we obtain 



i d 2 a 2 d 2 a\d 2 aa , _ . 



(ai ~ U{a 2 ,t 2 ,ai,t 1 ) 

xU(ai,ti,a ,t )p(ao,to), (81) 



cf. endnote [56| . Integration by parts was used to 
move the derivative to U(a 2 ,t 2 , ai,ti); square brack- 
ets emphasize that the differentiation does not apply to 
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U(ai, ti,ao,to). Similar considerations yield 

/ 7, s ~, . \ f d 2 otod 2 : Qi d 2 an 



x U(a 1 ,t 1 ,a ,t )p(a ,t n ), (82) 



(A\tx)S>(t2)) = 



d 2 OL 2 d 2 OL\d 2 Un 



1 d 

xU(ai,ti,a ,to)p(a ,to), (83) 



Y(a 2 ) 



(y(t 2 )AHh)) 



d 2 a 2 d 2 a\d 2 an 



Y{a 2 ) 



(at 



1 d \ 
+ 2~da~l) U ( a2,t2,ai,tl " > 

x U(at, t 1 ,a ,t )p(a Q ,t Q ). (84) 



We remind the reader that Eqs. ( [H])-([EI]) hold if t < 
t\ < t 2 . The latest operator in them is in fact arbitrary. 

Equations ([8Tj) - (|84p are expressions of "generalised 
phase-space correspondences" discussed in section III Bl 
They are exact and not associated with the path-integral 
representation of the phase-space amplitude. However 
their most natural interpretation is in terms of path- 
integral averages, with the derivatives related to "quan- 
tum jumps" of the trajectories. 



IV. ANALYTICAL EXAMPLE: THE KERR 
OSCILLATOR 

As a simple illustrative example we apply the "gen- 
eralised phase-space correspondences" (f2"U|) to the Kerr 
oscillator introduced in section [TTJ The same model was 
used in Ref. [24| to illustrate the effect of quantum correc- 
tions to the -W picture. This problem is exactly soluble; 
better still, all calculations implied by Eq. ( |2"0")) may be 
completed analytically. This makes the Kerr oscillator 
an ideal first testing ground for our approach. 



We assume that the oscillator is initially in a coher- 
ent state (this setupclosely mimics the collapse-revival 
experiment of Ref. 53[): 

<-> = </3| 



\(3), d\(3) = m- 



(85) 



The time-normally ordered correlation function is then 
easily calculated: 

GnihM) = {AHh)A{t 2 )) 

= |/3| 2 exp{|/3| 2 [ e - lK(t2 -* l) -l]}. (86) 

The subscript "H" distinguishes this as an exact Hilbert- 
space result. This expression follows from the exact so- 
lution for the Heisenberg operators, 



A(t) 



-t into* a 

a' e , 



(87) 



Gu(h,t 2 ) = (/3|aV* lKats e- 4 * 2Ka+a a| 



so that 



Equation (|86[) is found by expanding the exponents in 
power series and recalling the expansion of a coherent 
state over the number states. 

Calculations associated with Eq. ( |2"0")) take more ef- 
fort but are also quite straightforward. The equation 
for phase-space trajectories in the truncated Wigner rep- 
resentation is given by (fl"5]) with s(t) = 0. Phase- 
space evolution only affects the phase of a(t), so that 
\a{t)\ 2 = \a{0)\ 2 . With this observation Eq. ( [15]) is 
solved trivially, 



a(t) = a(0)exp{ - int(\a(0)\ 2 - l)}. 



(89) 



Stochasticity only enters the picture through the initial 
condition for a(t), distributed with probability 



W(a,a* 



■exp{-2|a(0)-/?| 2 }. 



(90) 



Strictly speaking, W(a,a*) is the Weyl-ordered 
quasiprobability distribution, or Wigner function, of the 
state |/3), but with W(a,a*) > such formal niceties 
may be disregarded. 

By making use of Eqs. ( I8"9"| and (f9"0"| we find for the 
symmetrically-ordered correlation function: 



Gw(*i,*a) = a*(ti)a(t2) 

For t x = h = we have G w (0, 0) = \(3\ 2 
leading to, 



A 2 + \- l -f (h-t 2 ) 

[i-f(tx-t 2 )Y 



■ exp 



in(ti - t 2 )- 



2 -i + f(h-t 2 ) 

1-f (tx-t 2 ) 



(91) 



1/2 as expected. The response terms in Eq. ( I2U1) are also easily calculated, 



da(t 2 ) 1 + f (t!-i 2 )(2|/3| 2 -l) 



da(tx) 



[l-f(ti-f 2 )]' 



• exp 



in(ti - t 2 )- 



2 -l + i±( tl -t 2 ) 
l-f(h-h) 



, t 2 > tx- 



(92) 
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FIG. 1: Normally-ordered correlation function for the Kerr-oscillator. Top and middle rows: Comparison between the exact 
solution Gn(ti,t2) (dashed line), the truncated- Wigner solution Gn(ti,ts) (solid line) and the "naively corrected" solution 
G c (ti,t 2 ) (dotted line) for |/3| 2 = 1, 2, 4, 8 (from left to right). Top row of graphs depicts the scaled moduli and the middle 
row — scaled phases, cf. Eqs. ( 196(1 . (|97[l . Bottom row: relative erros of various approximations to the normally-ordered 
correlation function (|98[) . Truncated-Wigner approximation Gn^i,^) (solid lines), uncorrected (symmetric) truncated- Wigner 
solution Gw(ti,t2) (dash-dotted lines), and the "naively corrected" solution Gc(ti,t2) (dotted lines). All quantities are plotted 
against the scaled time difference |/3|«;At. 



The quantity da(t\) / dafe) for t\ > t 2 is given by the same expression. Unlike Eq. ( [86)) which is exact, Eqs. ( |9T|) 

and (f9"2")l are approximations within the truncated Wigner approach. Combining them, we find the truncated- Wigner 
approximation to the normally ordered correlation function, 



[i-f (h-h)] 



■ exp 



in(ti — t 2 ) 



|/?| 2 -l + f ( tl -t 2 ) 



1 



K*i -*a) 



(93) 



Both the exact formula and the approximate formula for the time-normally ordered correlation function depend on 
the time difference At = t\ — 

We can demonstrate the accuracy of this approximate result by considering the series expansion, 



log [G N (At)/G H (At)] 



K 2 At 2 



12 



| 2 + 1) « 3 Ai 3 
I 



1 

9G 



3) K 4 At 4 + O (At 5 ) , 



(94) 



from which we see that Eq. (|93p is a good approximation 
if n\At\ <C 1, |/3|K|Ai| ~ 1. In other words, it holds over 
the collapse time scale, \At\ ~ 1/|/3|k, but fails over the 
revival time scale, \At\ ~ 1/n. 



In Fig. [T] we compare the exact function Gh(^i,^2) 
to its truncated- Wigner approximation Gn^i^) and to 
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the "naively" corrected correlation function 
G (ii,t 2 ) = GSv(*i,ia)-l/2. 



(95) 



(the latter corresponds to using the free-field commutator 
in place of the Heisenberg one) for different values of (3. 
Plots in the top and middle rows depict, respectively, the 
scaled modulus 



g(t 1 ,t 2 ) = \G(h,t 2 )\/\P\ 2 
and the scaled phase 

0(ti,ta) = argG(ti,*2)/|/?|7r 



(96) 



(97) 



of the correlation function 
Gn(ti,t2), Gn(*1) h), Gc(ti,t 2 ). 
torn row show the relative error, 



for G(*i,*a) 
Plots in the bot- 



S(h,t 2 ) 



G(ti,t 2 ) 



(98) 



Here we also include the symmetric correlation function 
G(ii, *2) = Gw(ti, t 2 )- Each column of plots corresponds 
to one value of /3: from left to right, \f3\ 2 = 1, 2, 4, and 8. 
All graphs are plotted against the scaled time difference 
\P\nAt. 

We see that the accuracy of Eq. (f9"3"]l is always supe- 
rior to that of the uncorrected as well as the naively 
corrected symmetric average. The response correction 
brings the truncated Wigner prediction into excellent 
agreement with the true solution for |/3| 2 > 2, but is 
not very accurate for \(3\ 2 = 1. The lack of accuracy for 
|/3 1 2 = 1 is not an unexpected result as the truncation 
process is generally thought to be justifiable as long as 
the number of quanta is significantly greater than the 
number of modes [54j . What is perhaps surprising here 
is how accurate the approximation becomes for |/3| 2 as 
small as 2, although we must remark that accuracy in 
calculating one particular operator moment does not im- 
ply accuracy in the calculation of all possible moments. 
It is worthy of reminding the reader that inaccuracies in 
Eqs. ( [9lJ (J93j) are solely due to the approximate nature 
of the truncated Wigner approach. By itself, Eq. ( [217)1 is 
exact. 



V. NUMERICAL EXAMPLE: THE 
BOSE-HUBBARD CHAIN 

In this section we apply our method to the one- 
dimensional Bose-Hubbard model. This model describes, 



in particular, neutral bosons in deep optical lattices and 
can be readily realized in experiments (see Ref. (55| for 
an overview). The Bose-Hubbard model is described by 
the Hamiltonian, 



N r 



H 



n J2 

k=l 



+ -n k (hk - l) 



- J(a\a k 



(99) 



with fik — a' k ak and fife, a k being the standard creation- 
annihilation pair for the k-th site, 



[ak, 



Skk' 



(100) 



with Skk' being the Kronecker delta. For simplicity we 
will consider a closed ring with periodic boundary con- 
ditions. The indices are to be understood modulo N, 
a N+i = Q-i- All definitions given in section III Al for the 



Kerr oscillator apply with the replacements 



k ■ 



A(t) -> A k {t), and A ] {t) -> A\{t). 



■ a k , 
The 



truncated- Wigner equations of motion for the system de- 
scribed by the Hamiltonian (|99|) read 



ictk 



c(|afe| 2 - l)ak - J(ak+i + a k -i) ■ (101) 



These equations follow from the traditional phase-space 
methods [20j and may be generalised to multitime aver- 
ages using the path-integral approach in section ITTT1 

Our aim is now to calculate the two-time normally or- 
dered correlation function, 



G 



Hfcfc 



(U,t 2 ) = (AlitJAk'ih)). 



(102) 



For two or more sites, this is a real problem: the prob- 
lem is not exactly soluble, nor can the calculations in 
the Wigner approach be done analytically. For the lat- 
ter, a natural choice is numerics in phase-space; after all, 
making real systems amenable to such methods is our 
ultimate goal. We employ the obvious generalisation of 
the one-mode formula (|20p : 



14 



a* k (h)a k > (t 2 ) 



1 da k > (t 2 ) 

2 0a fc (ti) 



a fc(*i)afc'(*2) 



da k ,(t 2 ) 



ti < t 2 , 



ti>t 2 . 



(103) 



We write this formula as an approximate one implying 
that the numerics are done with the truncated Wigner 
representation. In this case the averaging on the RHS 
reduces to that over the initial Wigner function, while the 
trajectories obey Eqs. ( UOip . Were the bar replaced by 
the double bar denoting a full path-integral quasiaverage, 
Eq. ( 1103)) would become exact. 

Importantly, implementing Eq. ( I103[) does not require 
independent "quantum jumps" at every time step. In 
fact a jump at zero time suffices. For each trajectory one 
can then use the chain formula 



da k < (h) 



E 



da k (t 2 ) da k »(t ) 
da k »{t ) datk'ih) 

da k (t 2 ) da* k „(to) 
da* k „{t ) da k '{h) 



(104) 



with to = 0. The quantities da k "(to)/da k '(ti) 
and da k „(to)/da k '(ti) are found by inverting the 
matrix comprising da k (h) / da k > (t ) , dal(ti) / da k > (t Q ) , 
dak(ti)/da%,(t ), and da>l(ti)/ da%, (t ) ■ Further de- 
tails can be worked out by complementing Eq. ( 
I104j) with similar chain relations for da k (t 2 )/da k '(ti), 
da k (t 2 )/da k ,(ti) and dal(t 2 ) / da k ,(ti), and using 



da k i(tx) 
da* k (h) 



da* k ,(h] 
da k (ti) 



= 8, 



kk' 



(105) 



da k >(h) da%,{ti) 



= 0. 



Numerical implementation of Eq. ( I103[) thus requires a 
minimum of 2N + 1 trajectories run in parallel, for every 
initial condition generated from the distribution 



JV N 



W(a(p),a*(p))= (-) n exp [ - 2|« fe (0) - /3 fc | 2 ]. 

(106) 



k=l 



This formula implies that the initial condition (the 
Heisenberg p-matrix) we use when evaluating (|103[) is 
a direct product of coherent states, 



|/3)= \l3 1 )®\f3 2 )®---®\(3 N ). 



(107) 



For better numerical performance we implemented four 
independent shifts per mode, requiring 47V + 1 trajecto- 
ries per "coin toss." Such numerical cost is obviously not 
prohibitive. 



As we wish to have a "reference point" against which 
to compare our results, we can use either exact diagonal- 
ization, which forces us to limit the number of sites in the 
Bose-Hubbard chain to two or three or to use the time- 
evolving block decimation algorithm (TEBD) 19]. The 
latter assumes small entanglement in the chain which is 
justified for not too large times and sufficiently small sys- 
tems. In Fig.[2]we plot the result of the truncated- Wigner 
calculation of the normally-ordered correlation function 
for N s — 10 sites and compare this to TEBD simulations. 
As the TEBD algorithm favors open boundary conditions 
we here (and only here) use these conditions. The figure 
shows the scaled modulus, 



ffNfcfe' (ti,t 2 ) 



|GNfcfc' *2)| 



(108) 



(left), and the scaled phase, 



, r , , x arg G Nkk >{t 1 ,t 2 ) . , 

<PNkk'{ti,h) = r-51 , (109) 



(right), of the correlation function for k = 5 and k' rang- 
ing from 1 to 10. In other words, each line in Fig. [^cor- 
responds to correlations between site 5 and either itself 
or some other site. All quantities are plotted versus the 
scaled time difference |/?|KAt. The inital condition was 
chosen as the same coherent state \(3 k ) in all modes with 
f3 k = \/2, k = 1, • • ■ ,10 (i.e., two quanta per mode). The 
hopping strength was set to J = 0.1 and the interaction 
strength to K = 1. The time t 2 was chosen arbitrarily as 
t 2 = 0.45. The average for the truncated- Wigner method 
was over 80, 000 runs. 

One recognizes a rather good agreement. Since the 
TEBD calculations are expensive we will resort in the 
following examples to the case of two and three modes, 
where direct numerical calculations in the full Hilbcrt 
space remain doable. 

In Fig. [3] we compare the results of the phase-space sim- 
ulations (solid lines) to those in the Hilbert space (dashed 
lines) for the two mode case. The plotted quantities 
are gn(h,t 2 ) (top left), gi 2 (ti, h) (top right), 0n(ti, t 2 ) 
(bottom left), and <fii 2 (ti,t 2 ) (bottom right). That is, 
the top row shows the modulus while the bottom row — 
the phase of the normally-ordered correlator; the left col- 
umn depicts the same-site, while the right column — the 
neighbour-to-neighbour correlations. Dependence of all 
quantities on t\ and t 2 is expressed naturally by 3D plots. 
However, while t\ changes continuously, t 2 is limited to 
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FIG. 2: The normally-ordered correlation function for the Bose-Hubbard chain of 10 sites simulated using the truncated- 
Wigner approximation with response correction (solid line) and the TEBD-method (dashed line) for comparision. Left: the 
scaled modulus g^^y (£i, fa), right: the scaled phase <^Nfcfc' [tit fa), cf. Eqs. ([TDSJ, (|109[1 . for k = 5 and k' = 1, • ■ ■ ,10. The 
initial condition (3k = k = 1, • • • ,10; J = 0.1 and k = 1. Graphs are plotted versus the scaled time |/3|kA£ with fa chosen 
arbitrarily as fa — 0.45. 



discrete values, t% = 0.1, 0.2, . . . , 1.8. The initial condi- 
tion is a coherent state \0) with = \[2 in each mode, 
the hopping strength is set to J = 0.1 and the interac- 
tion strength to k = 1. The average is over 80,000 runs. 
We see that the conditions imposed by number conserva- 
tion, gu(t,t) = 1, 4>u{t,t) = 0, are clearly met in Fig. [3] 
Furthermore, the agreement between Hilbert space and 
the truncated- Wigner approximation is very good. The 
modulus of the correlation functions is always well repro- 
duced. The phase appears to be not as well reproduced, 
but this impression is deceptive. In fact, error in phase 
increases when the modulus becomes small. Even in this 
case, the truncated- Wigner approach gives a reasonably 
good approximation to the phase of the correlation func- 
tions. 

A detailed comparison between the truncated- Wigner 
and Hilbert-space calculations for three sites may be seen 
in Figs. 3] and [5] The initial condition in Fig.[4]was chosen 
as three identical coherent states, 

\0) = \0) ® |/?> ® \0), (Fig.iJ (HO) 

with = V2~, while in Fig. [5] the coherent states differ in 
phases, 

\p) = \0)®\0e 2i */ 3 )®\0e ii */ 3 ), (Fig.© (Ill) 

with the same 0. Calculations were performed for 
k = 1 and J = 0.1,1,10. In both figures, the top 
row of graphs shows the quantities 311(^1,^2), 512(^1^2) 
and Sis (ii, £2) plotted versus the scaled time difference 
|/3|/t(ti — £2)1 with t% chosen arbitrarily as t% — 0.45. 
The grouping of graphs in Figs. 2] and [5] with respect 
to the values of J is self-explanatory. The results of 
truncated- Wigner calculations using Eq. ( I103P are shown 
as solid lines, the dashed lines represent the Hilbert- 
space results, and the dash-dotted lines — the uncorrected 



(symmetrically-ordered) correlation function. The mid- 
dle row of graphs depicts the corresponding scaled phases 
0n(ti)*2), 4>n{t\it2) and 013(^1,^2)- The bottom row 
represents the relative error of the truncated- Wigner sim- 
ulation compared to the Hilbert-space result, 



<W(*l'*2j 



GNfc/c' (^1,^2) 



Gnkk' (*i)*2) 



(112) 



Phase space averages were taken over 80, 000 runs. 

Once again, we see that the requirements imposed by 
number conservation are met in all results. For t\ = ti, 
9\\{t\M) = 1 and (j>n(h,t 2 ) = 0. 

The importance and accuracy of the response correc- 
tion manifests itself pretty impressively for J = 10 where 
the uncorrected phase-space solutions oscillate. This is 
not a numerical artifact because oscillations occur ir- 
respective of the time step of the numerical integra- 
tion. These oscillations cancel out with similar ocilla- 
tions in the response term leaving the correlation function 
smooth and in good agreement with the Hilbert-space re- 
sult. 

Looking at the error plots one can recognize that the 
method performs well over the collapse time scale. It is 
worthy of stressing that we look at the relative and not 
at the absolute error. The error tends to get large for 
At — > 1/|/3|k but actually this is only due to the fact that 
the absolute value of the correlation function is already 
very small. 

To develop a feeling for the statistical error, in Fig. 
[S] we compare the discrepancy between results obtained 
from two independent phase-space simulations with the 
relative error of the same pair of results compared to 
the exact (Hilbert-space) result, for the two mode Bose- 
Hubbard chain. The initial condition is a coherent state 
with — y2 in each mode. The hopping strength is 
set to J = 0.1, and the interaction strength to k = 1. 
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FIG. 3: Comparison of the truncated- Wigner results (solid lines) to Hilbert-space ones (dashed lines) for the normally-ordered 
correlation function of the Bose-Hubbard chain of two sites. Top row: the scaled moduli (?n (fi, £2) (left) and <?i2(ii,i2) (right), 
bottom row: the scaled phases 0n(ti,ta) (left) and ^i2(ii,ta) (right). The initial condition is /3i = fi% = J = 0.1 and 
k = 1. The time t\ changes continuously while t% is limited to discrete values ranging from 0.1 to 1.8. 



The relative error between two correlation functions sim- 
ulated in phase-space with the same number of "coin 
tosses" is shown as solid lines. The dashed lines are used 
for errors of the same pair of simulations in phase-space 
compared to the Hilbert-space solutions. The picture on 
top shows the said errors for the on-site average with 
k = k' = 1, the one at the bottom for the averages be- 
tween two modes, k = 1, k! = 2. All graphs are plot- 
ted versus the scaled time with ti arbitrarily chosen as 
t 2 = 1.3. 

It is evident from the figure that the statistical er- 
rors are insignificant. The discrepancy between results 
in phase-space are either just small, or small compared 
to the accuracy of the method. At times \/3\nAt where 
the relative errors of the phase-space result compared to 
the Hilbert-space solution are of the same order as the 
errors of two independent runs in phase-space, the corre- 
lation function is already very small. 



All in all, the "response correction" brings the results 
of the truncated- Wigner simulation in phase-space into a 
good agreement with the Hilbert-space simulations. Our 
observation for a single mode, that is, that the method 
gives good results over collapse time scales but fails on 
revival times, also holds for two and three modes. It is 
worthy of reminding the reader that all errors are solely 
due to the approximate nature of the truncated Wigner 
simulation. By itself, Eq. ( 1 1 03f> is exact. 



VI. SUMMARY 

A phase-space path-integral approach generalising the 
symmetric representation of Schrodinger operators to the 
Heisenberg picture is developed, and "generalised phase- 
space correspondences" allowing one to commute Heisen- 




|/3|«At 

FIG. 4: The normally-ordered correlation function for the three mode Bose-Hubbard chain for the initial condition (|110[) . k = 1 and J = 0.1, 1, 10. Top row: the 
quantities <?ii(£i, £2), 312(^1^2) and gri3(ii,i2); grouping of graphs with respect to the values of J is self-explanatory. Truncated- Wigner calculations (solid lines), 
Hilbert-space results (dashed lines), and the uncorrected (symmetric) correlation function (dash-dotted lines). Middle row: the corresponding scaled phases </>ii(ii, £2), 
</>i2(ii,t2) and 0i3 (tx, t<z). Bottom row: the relative error of the truncated- Wigner simulation compared to the Hilbert-space result. All quantities are plotted versus 
the scaled time difference |/3|fi;Ai with t% choosen arbitrarily as ti — 0.45. 
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berg operators with unequal time arguments are de- 
rived. The conventional truncated Wigner representation 
emerges as an approximation within the path-integral ap- 
proach. This results in formal techniques allowing one to 
calculate time-normal averages of Heisenberg operators 
approximately with relative ease. These techniques have 
been verified for the Kerr oscillator and for the Bose- 
Hubbard model showing a good agreement with exact 
Hilbert space calculations at collapse time scales for sur- 
prisingly low numbers of oscillator quanta. 
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FIG. 6: Assessment of statistical errors for two mode case 
with J = 0.1, k = 1, \0) = 1/2(1, 1). Solid lines show the dis- 
crepancy between results derived from two independent runs 
with the same number of "coin tosses" (40,000 and 80,000). 
Dashed lines show the relative error of the same results com- 
pared to the exact (Hilbert-space) results. All graphs are 
plotted versus the scaled time \(3\nAt, with ti set arbitrarily 
to ti — 1.3. We see that the statistical errors are insignificant. 
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